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Abstract. The overlapping stage of heavy-ion reactions can be simulated by dynamical mi- 
croscopical models, such as those built on the basis of the Molecular Dynamics (MD) ap- 
proaches, allowing to study the fragment formation process. The present performances of the 
Quantum MD (QMD) code developed at the University of Milano are discussed, showing 
results concerning fragment and particle production at bombarding energies up to < 700 
MeV/A, as well as a preliminary analysis on the isoscaling behaviour of isotopic yield ratios 
for reactions with isospin composition N/Z in the (1 - 1.2) range, at a 45 MeV/A bombarding 
energy. 



1 Introduction 

In the phase diagram of nuclear matter two different regions, where signals related 
to the probable onset of phase transitions have been detected, have been explored 
with increasing interest. First of all, a phase-transition is expected to occur at a 
critical temperature of a few MeV and a critical density relatively low compared 
to the normal nuclear matter density. Subsaturation densities can be accessed by 
heavy-ion collisions in the intermediate energy domain (30 MeV/A - 200 MeV/ A), 
which allow to study the thermodynamic properties of strongly interacting nuclear 
matter. The characteristic features of nuclear forces (long-range attraction and short- 
range repulsion) could be responsible of the so called liquid-gas phase transition in 
this region. Symmetry energy effects and Coulomb effects are recognized to play an 
important role, whereas the role of the nucleon spin has still to be investigated [1]. 
Reactions in this energy domain are currently studied by means of a number of 
facilities worldwide. In particular, we mention the experiments performed by the 
INDRA Collaboration at Ganil, by the ALADIN Collaboration at GSI, and by the 
CHIMERA and ISOSPIN Collaborations at the LNS in Catania. On the other hand, 
relativistic and ultrarelativistic energy heavy-ion collisions can be used to explore a 
region at much higher (p, T) where another phase transition probably occurs: the one 
from hadronic matter made of baryons and mesons to the quark-gluon phase, already 
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investigated at RHIC, and which will be further explored by the LHC heavy-ion 
program. According to lattice calculations, the critical temperature for the last phase 
transition or sudden cross-over into a deconfined quark-gluon plasma, is expected 
to be « 180 MeV. In this paper we will focus on the first kind of reactions, with 
projectile energies in the range from tens MeV/A to a few hundreds MeV/A. 

Among the possible signatures of a liquid-gas phase transition, we mention in 
particular multifragmentation, i.e. the simultaneous breakup of an excited nuclear 
system in a large number of IMF (Intermediate Mass Fragments, 3 < Z < Z max , 
where Z max is well below the total charge of the nuclear system) (liquid phase), 
coexisting in case with LMF (Light Mass Fragments, Z < 3) (gas phase). Multifrag- 
mentation is a universal phenomenon observed in intermediate-energy nuclear in- 
teractions induced by hadrons/photons/heavy-ions [2] . During the overlapping stage 
of heavy-ion collisions (typical time w 100 - 200 fm/c), depending of the projectile 
ion bombarding energy and impact parameter, the nuclear system can undergo com- 
pression and reach high excitation energies. As a consequence, it starts to expand 
and can go on expanding down to sub-saturation densities (p w 0.1 — 0.3 po, where 
po is the normal nuclear density) and reach temperatures T « 3 — 8 MeV, where 
it becomes unstable and breaks up into multiple fragments. These conditions are 
typical of the liquid-gas mixed phase [3]. 

One of the issues still open is if a thermal equilibration is reached in these re- 
actions. The statistical description of multifragmentation relies on this assumption. 
According to this approach, the emissions from multifragmenting sources are not af- 
fected by the dynamics which have led to source formation. Only the source global 
thermodynamical properties are important to determine the fragment pattern. The 
dynamical description of multifragmentation instead does not rely on this assump- 
tion. According to the dynamical approach, multifragmentation is a fast process, and 
the IMF detected as reaction products are related to nucleon-nucleon correlations 
already occurring during the ion-ion overlapping stage and surviving the following 
expansion phase (memory effects) [4]. 

The fast stage of heavy-ion collisions can be simulated by means of Monte Carlo 
dynamical models, such as the Quantum Molecular Dynamics ones (QMDs), allo- 
wing to describe phase-space fluctuations in the ion-ion overlapping stage, which 
lead to fragment formation [5]. These are microscopical models, at nucleon level. 
Preferring a microscopical model to a macroscopical one is mainly motivated by the 
fact that a large number of degrees of freedom can appear in a reaction. Depending 
on the bombarding energy of the projectile system and on the impact parameter, 
one has to consider that excitation and deformation of projectile and target ions, 
neck formation and breaking, nucleon transfer, different way of fragmenting the 
composite system, and many other effects can simultaneously play a role [6]. It 
could be very difficult to take into account so many degrees of freedom in a macro- 
scopical model. In a microscopical dynamical model, on the other hand, they can be 
described in a more natural, automatic and unified way. 

At the end of the fast stage ( w 10~ 22 s) of ion-ion interactions simulated by 
QMDs, excited fragments can be present. In order to reproduce experimental data 
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concerning particle and fragment emission, also the de-excitation process has to 
be accounted for. Fragment de-excitation takes place on a time scale (« 10~ 21 - 
10~ 15 s) far larger than the ion-ion overlapping process and is thus more precisely 
and faster described by other models. Anyway, some attempts to extend dynamical 
simulations up to a few thousands fm/c (i.e. < 10 -20 s) have also been carried out 
with success by means of modified advanced versions of QMD (see e.g. [7]). 

In this paper, the QMD code developed in Milano [8, 9] is used to simulate 
the overlapping stage of ion-ion collisions, leading to excited fragments (prefrag- 
ments). Primary fragments are identified at the end of each QMD event (t « 200 
- 250 fm/c) by applying a phase-space correlation algorithm to the nucleon distri- 
bution. Two particles are supposed to belong to the same cluster if their relative 
distance and momentum are within fixed amounts. Our QMD code has been inter- 
faced to the fragment de-excitation model already implemented into the general pur- 
pose FLUKA transport and interaction code [10-12]. The evaporation/fission/Fermi 
break-up module [13] available in FLUKA allows to describe primary fragment de- 
excitation on a statistical basis. General features of both codes have been presented 
in previous papers. In this paper, results of simulation concerning particle and frag- 
ment emission at energies within a few hundreds MeV/A are shown and compared 
to experimental data available in literature. In particular, the code has been applied 
both to the simulation of heavy-ion collisions up to < 700 MeV/A, as described in 
Section|2]and[3] and to the study of isospin effects at energies w 50 MeV/A, typical 
of multifragmentation events, as explained in Section|4] Our conclusions, with some 
reference to practical applications too, are given in Section|5] 

2 Particle production 

As far as thin targets are concerned, double-differential neutron production cross- 
sections for inclusive reactions at energies around a few hundreds MeV/A have been 
presented in Ref. [8]. Target thickness has been neglected in our simulations. 

As far as thick targets are concerned, a recent paper [14] claimed that old data on 
neutron double-differential yields, published in previous papers [15, 16] by the same 
working group, suffer from an underestimation of the neutron detector efficiencies, 
with errors increasing with increasing neutron energy. The authors of [14] have thus 
re-evaluated the old results taking into account updated estimates of the detection 
efficiencies and published new data from recent measurements. We simulate the 
same reactions, taking into account target thicknesses, chosen in the experiment 
in such a way that the energy losses of projectile ions lead to their stop inside the 
targets. The results of our theoretical simulations of neutron emission from C and Ne 
beams impinging on C, Al and Cu systems at a 400 MeV/A bombarding energy are 
shown in Fig. [T] together with both the re-evaluated and the old experimental data, 
where available. It is apparent that, even if, in general, the results of the theoretical 
simulations better agree with the re-evaluated data than with the old ones as far 
as the tails at higher energies are concerned, there are energies/angles where the 
theoretical simulations better agree with the old data than with the re-evaluated ones. 
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In particular, in case of C projectiles, the results of the simulations always overe- 
stimate the experimental data for the emission of high-energy neutrons at a ^ 60° 
angle. The same trend at 60° has been retrieved by using a completely independent 
code, e.g. the PHITS code [18] including the JQMD model [19], as shown in [14]. 
These facts point out that, besides the uncertainties indeed present in the models for 
theoretical simulations, further uncertainties in the experimental data could exist. 
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Figure 1. Double-differential neutron production yield for C (left) and Ne (right) projectiles 
impinging on C (top), Al (center) and Cu (bottom) at a 400 MeV/A bombarding energy. The 
results of the theoretical simulations made by QMD + FLUKA are shown by solid triangles, 
whereas the experimental data re-evaluated by Satoh et al. [14] are shown by filled circles. 
The old experimental data [17] are shown by empty circles. In each panel, distributions at 0°, 
7.5°, 15°, 30°, 60°, 90° angles with respect to the incoming beam direction in the laboratory 
frame, have been multiplied by decreasing powers of 10, for display purposes. 
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3 Fragment production 
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Figure 2. Projectile-like fragment production cross sections for 650 MeV/A Ar beams frag- 
menting on C (top left), Al (top right) and Ar (bottom) targets, calculated by QMD + FLUKA 
(solid histograms) in comparison with experimental data taken from Ref. [20] (points). For 
Z = 1 only the contributions to the cross-sections due to fragments (A > 1) are shown. 



As far as fragment production is concerned, we calculate projectile-like frag- 
ment production cross-sections for CI and Ar projectiles impinging on C, Al, Cu 
and Pb targets at a 650 MeV/A bombarding energy. Experimental data concerning 
these reactions were taken at the NSRL at the Brookeven National Laboratory, US 
and at the HIMAC at the National Institute of Radiological Science, Japan and pu- 
blished in Ref. [20]. The results of our simulations are shown together with the 
experimental data for Ar and CI projectiles in Fig.|2]and [3] respectively. Despite the 
overall reasonable agreement, it is apparent that the theoretical simulations under- 
estimate the production of the most massive fragments, in case of the most massive 
targets. These fragments are emitted in very peripheral collisions. A reason of this 
underestimation can be ascribed to electromagnetic dissociation, which increases 
with the system mass and was not accounted for by our simulations. Anyway, we 
estimated that this effect is not enough to explain the discrepancies. On the other 
hand, for lighter fragments, corresponding to less peripheral collisions, the agree- 
ment between theoretical expectations and experimental data is more reasonable and 
the odd-even effect observed in experimental data is reproduced by the simulation. 
In our simulations the last feature can be ascribed to de-excitation emission effects. 
For the lightest fragments which can be identified by the experiment the agreement 
is sometimes slightly worser. The experimental procedure allows charge identifica- 
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tion by means of the AE deposited in silicon detectors. With decreasing charge, 
the peaks in the (Z, count of events) plots which allow charge identification become 
increasingly broader [20]. For Zf rag < Z pr i mary /2 it was not possible to distinguish 
the peaks corresponding to different charges with the large-acceptance (~ 5° — 13°) 
detectors used in the experiment. Thus, experimental data were obtained for frag- 
ments in the range 5 < Zf rag < Z pr i mary /2 by using small acceptance (~ 1° — 2°) 
detectors. Furthermore, while the detectors measure fragmentation within a small 
angle around the beam axis, the angular distributions of lighter fragments are less 
peaked in forward direction. Thus, in our opinion, the uncertainties in charge identi- 
fication together with those on the estimation of the correction due to the broadness 
of the angular distributions could largely affect the results on the cross-sections for 
light fragment production. 
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Figure 3. Projectile-like fragment production cross sections for 650 MeV/A CI beams frag- 
menting on C (top left), Al (top right), Ar (bottom left) and Pb (bottom right) targets calcu- 
lated by QMD + FLUKA (solid histograms) in comparison with experimental data taken from 
Ref. [20] (points). For Z = 1 only the contributions to the cross-sections due to fragments (A 
> 1) are shown. 

The same trend has also been observed in case of different projectiles. The re- 
sults for Fe + H, Fe + C and Fe + Al at 400 and 500 MeV/A bombarding energies 
are plotted in the three panels of Fig. |4] The agreement between the results of our 
theoretical simulations and the experimental data [21] is well within 30 - 35 % for 
most fragments. The QMD + FLUKA interface gives more reliable results in case 
of heavy targets than in case of the hydrogen one. In case of the H target, the largest 
discrepancies have been obtained for the lightest fragments (Z = 12 - 16) detected 
by the experiment (the threshold for charge identification in these experiments was 
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Figure 4. Scatterplots of fragment production cross-sections predicted by QMD + FLUKA 
de-excitation (vertical axis) vs. experimental cross-sections [21] (horizontal axis) for Fe ions 
impinging on H (top left), C (top right) and Al (bottom) ions at 400 (crosses) and 500 
(squares) MeV/A bombarding energies. The errorbars on the theoretical results are statistical 
only. The two lines respectively above and below the diagonal correspond to deviations of 
the modelled results from the experimental ones by ± 33%. 



Z = 1 2) at the lowest bombarding energy, whereas a better agreement between theory 
and experiment is observed for the Z = 22 - 25 fragment production cross-sections. 
In case of C and Al targets, the best agreement is observed for the Z = 14-21 
fragments, whereas the largest discrepancies are seen for the Z = 24 fragment cross- 
section, overpredicted by the theory. The question about the dependence of these 
results on the target initial configuration need further investigation. Up to now, we 
use only a few selected projectile and target ion configurations to perform our simu- 
lations, giving rise to different events by rotating in a random way the coordinates of 
these nuclei. The simulations can be carried out as well by choosing a larger number 
of initial configurations (always satisfying some criteria concerning mean square ra- 
dius, binding energies, etc.) and one can study if this affects the results obtained for 
the most peripheral collisions. In fact, the most peripheral collisions are expected to 
be more sensitive than the central ones to properties of the initial ion configurations 
such as root mean square radius oscillations. 

Further effects are still lacking in our and in other QMD models. Among the 
others, we mention shell effects, a spin-orbit term in the Hamiltonian, ground state 
deformation of initial nuclei. In particular, accounting for shape deformation can be 
important for very heavy systems, since fusion/fission reactions are sensitive to the 



8 



M.V. Garzelli 



geometry of the collision and the relative orientation of the deformed nuclei. While 
the last effects can indeed be included by means of a proper choice of the initial 
nuclear states, introducing shell effects into a transport model is still a task of great 
difficulty [6]. 

Furthermore, a proper treatment of antisymmetrization is still lacking in our 
model. However, it has already been included in some antisymmetrized advanced 
versions of QMD like AMD [22] and FMD. These ones are slower, and have mainly 
been used to study nuclear structure properties, especially for light nuclei, even if 
several examples of simulation of nuclear reactions with AMD have recently been 
produced. 

4 Isoscaling observables 

The isoscaling technique allows to obtain useful information concerning the symme- 
try energy of isospin asymmetric nuclear matter, one of the quantities most uncertain 
and most discussed. This information is important to understand the mechanisms 
which regulate neutron star formation and structure, collapses of massive stars and 
supernova explosions. Both in theoretical simulations of nuclear collisions and in 
experiments the isotopic yield ratio from two reactions with similar total sizes and 
temperatures, but different isospin asymmetries (N — Z)/(N + Z), has been found 
to follow an exponential law [23,24]: 



where the index 2 refers to the reaction with the largest neutron number (more asym- 
metric) and the index 1 refers to the reaction with the largest proton number (more 
symmetric). This relation is valid at fixed Z, limiting to the isotopes (N, Z) for which 
the isotopic yield distributions can be approximated by gaussians. This means that 
the logarithm of the yield ratio is linear in N, at fixed Z, with a slope given by 
a(Z). An analogous relation is valid at fixed N, meaning that the yield ratio is also 
linear in Z. This result is known as the isoscaling phenomenon and has been exten- 
sively studied [25,26]. Among the others, systematic experiments were performed 
at the Cyclotron Institute of TAMU and at the NSCL of MSU, analyzing the systems 
Ca + Ni, Ar + Ni, Ar + Fe at 25, 33, 45 and 53 MeV/A bombarding energies [27]. 
The experimental isotopic yield ratios for light fragments show a linear dependence 
on N at fixed Z in the logarithmic plane, with a(Z) well approximated by a constant 
a for each fixed bombarding energy and couple (1, 2) of reactions. The results have 
been interpreted by the authors of [27] both in the framework of the AMD dynami- 
cal model and in the framework of a statistical model. We have made simulation of 
the Ca + Ni, Ar + Fe and Ar + Ni reactions at 45 MeV/A using the QMD model 
interfaced to the FLUKA de-excitation module. The isospin compositions N/Z for 
Ca + Ni, Ar + Ni and Ar + Fe amount to « 1.04, 1.13 and 1.18, respectively. Our 
results for the isotopic yield ratios at the end of the whole simulation are plotted 



y( 2) (N,z) 

*a)(N,z) 



= cxp(a(Z)N + X(Z)) 
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in Fig. [5] In each panel, also the results for the lightest fragments (Z = 1, Z = 2) 
are included. It is apparent from our simulation that the linear dependence on N of 
the fragment yield ratios at fixed Z is a good approximation of the results in the 
logarithmic plane. Slight deviations from this linear dependence can be due to low 
statistics (« 10000 events distributed over all possible impact parameters were si- 
mulated for each reaction). By comparing the left and the right panel it is evident 
that larger slopes are found for the yield ratios from the couple of reactions with 
larger difference in the isospin compositions (N2/Z2 — N1/Z1). This behaviour is 
expected and in qualitative agreement with the experimental observations. However, 
we observe values of a which differ with Z. We calculate an average a value using 
the a(Z) values in the Z = 3 - 7 range. We obtain average a values rj 0.18 and rs 
0.31, for the Ar + Ni / Ca + Ni and Ar + Fe / Ca + Ni ratios, respectively. These 
values are larger that the ones observed experimentally. A reason of these discre- 
pancies is the fact that only the fragments detected around a 44° angle were selected 
for performing the experimental isoscaling analysis, whereas we consider all frag- 
ments, independent of the emission angles. This means that in the experiment the 
largest contribution to the selected fragments comes probably from central colli- 
sions, and these fragments are expected to include nucleons from both the projectile 
ion and the target one, whereas in our inclusive simulations IMF from projectile- 
like or target-like residues de-excitation are included. A preliminary analysis on 
our simulated events show that the a value is indeed affected by the choice of the 
impact parameter, and decreases significantly when selecting only the most central 
events. This behaviour is in qualitative agreement with Ref. [23], which claims that 
isoscaling is observed for a variety of reaction mechanisms, from multifragmenta- 
tion to evaporation and DIS events, with different slopes in the logarithmic plane. 
Furthermore, we can observe that the results of our analysis are quite sensitive to 
the number of isotopes included in the linear fit at fixed Z or, in other words, to the 
goodness of the gaussian approximation to the fragment isotopic distributions. The 
last ones are plotted for fragments with Z = 3 - 6 in Fig. [6] In each panel the isotopic 
distributions from the three considered reactions are plotted together. The isotopic 
distributions for the emission of light fragments turn out to show memory of the ini- 
tial isospin content of the system of reacting ions. When comparing the results for 
the three reactions, it is apparent that the reaction with the smallest asymmetry (or 
the smallest N/Z value) gives rise to isotopic distributions with larger tails at smaller 
N, whereas the one with the largest asymmetry (or the largest N/Z value) gives rise 
to larger tails at larger N. This result is in qualitative agreement with experimen- 
tal observations. This behaviour is expected to become even more evident at lower 
bombarding energies. 

When considering the free particles emitted in the fast stage of the collision, 
our QMD simulation suggests that the average yield of emitted protons is larger 
than the one of emitted neutrons in the Ca + Ni reaction, which has the lowest N/Z 
value. The other reactions exhibit the opposite trend for the collisions at low impact 
parameters. This behaviour is due to Coulomb force. In particular, according to our 
simulations, at pre-equilibrium, the yield of protons Y p emitted in Ca + Ni central 
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Figure 5. Fragment isotopic yield ratios obtained at the end of our simulations by QMD + 
FLUKA de-excitation, as a function of N for a 45 MeV/A beam energy. The ratios for the 
Ar + Ni / Ca + Ni and Ar + Fe / Ca + Ni couple of reactions are shown in the left and in 
the right panel, respectively. The different symbols correspond to different values of Z in the 
Z = 1 - 8 (from the left to the right of each panel) range. 
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Figure 6. Fragment isotopic distributions for Li, Be, B and C produced by Ca + Ni (stars), 
Ar + Ni (circles) and Ar + Fe (squares) reactions at a 45 MeV/A bombarding energy, at the 
end of our simulations by QMD + FLUKA de-excitation. 



collisions is a fraction w 20% larger than the yield of emitted neutrons Y n , whereas 
Y p is k 10% and 15% smaller than Y n for Ar + Ni and Ar + Fe central collisions, 
respectively. 

We have also investigated up to which extent the values of a are sensitive to 
the use of the final fragments (from the full simulation) or the use of the primary 
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excited fragments in performing the yield ratio analysis. We have found that, for the 
reactions considered in this study, the average value of a at the end of the overlap- 
ping stage of the collisions described by QMD, denoted by cthot, can be larger than 
the value of a at the end of the full simulation including de-excitation effects by no 
more than w 20%. 

Finally, we calculate the fragment asymmetry (Z/A)i; q 2 of the liquid phase for 
the three reactions at t = 250 fm/c. According to Ref. [23], the liquid phase is defined 
as composed by all fragments with A > 4. For central collisions in each of the three 
systems the average (Z/ A)i; q 2 at the end of the pre-equilibrium stage turns out to be 
lower then the corresponding value at t = 0, as already found for the same reactions 
at a lower energy by [27], with a difference between the values at different time 
more evident for the most symmetric reaction (Ca + Ni). The dependence of the 
average (Z/A)i; q 2 on time and impact parameter is under further investigation. 

5 Conclusions and practical applications 

A few examples of the application of the QMD code developed in Milano and inter- 
faced to the FLUKA de-excitation module have been produced, in the simulation of 
heavy-ion collisions both at a low bombarding energy (E 50 MeV/A) and at a few 
hundreds MeV/A beam energies, showing results which are in agreement, at least 
from a qualitative point of view, with experimental data and/or previous studies. The 
studies at low bombarding energies are mainly driven by theoretical issues and are 
aimed at gaining a better understanding of the role of isospin in multifragmenta- 
tion processes and in astrophysical mechanisms crucial to determine the evolution 
of supernovae and neutron stars, and can give insights on the optimization of the 
production of rare isotopes far from stability. On the other hand, the possibility to 
extend simulations by the same model even at higher bombarding energies allows 
to cover application purposes too. 

Among the practical application of QMD models we mention hadrontherapy, 
civil aviation and space radiation protection, and single-event effects in microelec- 
tronics. As far as hadrontherapy is concerned, ion beams (in particular C beams) are 
used to produce Complex Lesions, i.e. multiple breaks of the DNA double strands 
contained in cells. This effect can be exploited in order to kill tumours. This ap- 
plication requires the integration of nuclear models with radiobiological models. In 
particular, the ion biological effectivenesses turns out to be larger than the X-ray 
and proton ones. A few facilities built on the basis of these principles already exist 
(e.g. in Japan and in Germany, at GSI and in Heidelberg) or are under construction 
(e.g. the CNAO in Pavia, Italy). At higher energies, crews in space missions and in 
air-flight are exposed to Galactic Cosmic Rays (GCR). Spacecraft and aeroplanes 
have to be carefully designed in order to minimize GCR interaction effects [28]. 
Furthermore, functional upset of microelectronic memory devices (silicon chips) 
in space missions can also occur and deserves particular attention [29]. This phe- 
nomen can even take place at sea level, due to atmospheric neutrons (produced in 
the interactions of Cosmic Rays with air) in the 50 - 1000 MeV energy range. 
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